Numerical Computation Workshop

Nina Dokeva
Paritosh Mokashi
Rob Knapp
Oliver Reubenkonig

Wolfram Research Numerical Computation Development Team

Stretching

Make a histogram of the distribution function for |λ|/n where λ is the largest eigenvalue of nxn with uniformly distributed real entries for some different values of m.

Here is the case n = 1:

Histogram[RandomReal[1,10^5],100,"PDF"]

Simple way

distribution[n_, samples_, bins_] := 
Histogram[Table[Abs[First[Eigenvalues[RandomReal[1,{n,n}],1]]]/n,{samples}],bins,"PDF"]

s = 10^5;
b = 100;
Do[
Print[Apply[StringForm["n = `2`, timing = `1`",Round[10 #1]/10.,#2]&,AbsoluteTiming[plot=distribution[2^k,s,b];2^k]]];
Print[plot],
{k,1,5}]

Model

s = 10^5;
mdata = Table[
{n = 2^k,EstimatedDistribution[Table[Abs[First[Eigenvalues[RandomReal[1,{n,n}],1]]]/n,{s}],NormalDistribution[μ,σ]]},{k,1,5}]

sdata = Apply[{#1,#2[[2]]}&,mdata,{1}]

ListPlot[Log[sdata]]

Fit[Log[sdata],{1,x},x]

smodel[nm_] = Exp[% /. x->0] /nm

Table[smodel[2^k],{k,1,10}]

Warm up

To start up your numerical computation thinking, lets start with an exercise.

Given a set of points in 30 dimensions:

SeedRandom[59715];
data = RandomReal[1,{10^6, 30}];

and a test point:

tpoint = RandomReal[1,30];

Write a simple program that will find the element of data nearest to the test point

Nearest[data->Automatic,tpoint]

and then try to make the program run as fast as possible.

Attempt 1

myNearest1[data_,point_] :=Module[{dist,mindist=$MaxMachineNumber,index},
Do[
dist = Norm[point - data[[i]]];If[dist<mindist,mindist=dist;index=i],
{i,Length[data]}];
index]

myNearest1[data,tpoint]

Something I find very useful to test timings is the following function

SetAttributes[AverageTiming, HoldAll];AverageTiming[e_] := Module[{t = First[AbsoluteTiming[e;]],tries = 1},
While[t<1,
If[t < 0.001,tries *= 1000, tries = Ceiling[5*tries/(4 t)]];
t = First[AbsoluteTiming[Do[e,{tries}]]];
];
t/tries
]

AverageTiming[myNearest1[data,tpoint]]

AverageTiming[Nearest[data->Automatic,tpoint]]

With the kernel code under development now, we get

AverageTiming[Nearest[data->Automatic,tpoint]]

Attempt 2

We don’t actually need to use Norm -- using the square distance is equivalent and a tiny bit faster

myNearest2[data_,point_] :=Module[{diff,dist,mindist=$MaxMachineNumber,index},
Do[
diff = point - data[[i]];dist =diff.diff;If[dist<mindist,mindist=dist;index=i],
{i,Length[data]}];
index]

AverageTiming[myNearest2[data,tpoint]]

Attempt 3

Using vector arithmetic is potentially faster (but does require more memory)

myNearest3[data_,point_] :=Module[{diff,dist,mindist=$MaxMachineNumber,index},
diff = Map[point-#&,data];
dist = MapThread[Dot,{diff,diff}];
Ordering[dist,1]
]

myNearest3[data,tpoint]

AverageTiming[myNearest3[data,tpoint]]

Attempt 4

Try compiling:

myNearest4 = Compile[{{data,_Real,2},{point,_Real,1}},Module[{diff,dist,mindist=$MaxMachineNumber,index},
diff = Map[point-#&,data];
dist = MapThread[Dot,{diff,diff}];
Ordering[dist,1]
]]

myNearest4[data,tpoint]

AverageTiming[myNearest4[data,tpoint]]

Attempt 5

Sometimes the compiler does even better with simpler commands -- lets try compiling Do for saving memory:

myNearest5 = Compile[{{data,_Real,2},{point,_Real,1}},Module[{diff,dist,mindist=$MaxMachineNumber,index},
Do[
diff = point - data[[i]];dist =diff.diff;If[dist<mindist,mindist=dist;index=i],
{i,Length[data]}];
index]]

myNearest5 = Compile[{{data,_Real,2},{point,_Real,1}},Module[{diff,dist,mindist=$MaxMachineNumber,index=0},
Do[
diff = point - data[[i]];dist =diff.diff;If[dist<mindist,mindist=dist;index=i],
{i,Length[data]}];
index]]

AverageTiming[myNearest5[data,tpoint]]

Attempt 6

Sometimes the compiler does even better with simpler commands -- lets try compiling Do for saving memory:

myNearest6 = Compile[{{data,_Real,2},{point,_Real,1}},Module[{diff,dist,mindist=$MaxMachineNumber,index=0},
Do[
diff = point - data[[i]];dist =diff.diff;If[dist<mindist,mindist=dist;index=i],
{i,Length[data]}];
index],CompilationTarget->"C"]

AverageTiming[myNearest6[data,tpoint]]

The moral(s) of the story:

Using vector (listable) arithmetic can save lots of time
Using Compile with simple commands can be even faster
If you are going to use the function several times, using CompilationTarget->”C” can make it even faster yet.
and ...
If there is a Mathematica function to do the task, it may be much faster.

myNearest6 = Compile[{{data,_Real,2},{point,_Real,1}},Module[{diff,dist,mindist=$MaxMachineNumber,index=0},
Do[
diff = point - data[[i]];dist =diff.diff;If[dist<mindist,mindist=dist;index=i],
{i,Length[data]}];
index],CompilationTarget->"C"]

AverageTiming[myNearest6[data,tpoint]]

Aerobic set

Now to get warmed up to bigger things, another gentle exercise:

Write a simple univariate Newton root finder that will find the root of a function to a specified precision given a starting point.  (Assume the function can be differentiated symboilcally and that there is quadratic convergence at the root)

myRoot[f_,{x_,x0_}, p_] := ...

Some good functions to try it on are

f1 = x^2-2;x01 = 1; r1 = Sqrt[2];

f2 = Exp[x] - Sin[x]; x02=-4;Plot[f2,{x,-4,4}]

If you get something working and want to be more ambitious, write it to work so that the computation precision increases appropriately as the root is approached.

Attempt 1

print = Print;
maxcount = 10;

myRoot1[f_,{x_,x0_},p_] := Module[{df=D[f,x],r = N[x0,p],fx,dfx,count = 0,step,tol = 10^-p},
While[(count++)<maxcount,
fx = N[f /. x->r,p];
dfx=N[df /. x->r,p];
step = fx/dfx;
r = r - step;
print[{count,step}];
If[Abs[step] < Abs[r] tol,Break[]];
];
r
]

myRoot1[f1,{x,x01},30]

Attempt 2

myRoot2[f_,{x_,x0_},p_] := Module[{df=D[f,x],r = N[x0,p],fx,dfx,count = 0,step,tol = 10^-p},
Block[{$MinPrecision = p,$MaxPrecision=p},
While[(count++)<maxcount,
fx = f /. x->r;
dfx=df /. x->r;
step = fx/dfx;
r = r - step;
print[{count,step}];
If[Abs[step] < Abs[r] tol,Break[]];
];
r
]]

myRoot2[f1,{x,x01},30]

%-Sqrt[2]

Attempt 3

myRoot3[f_,{x_,x0_},p_] := Module[{df=D[f,x],r = N[x0],fx,dfx,count = 0,step,tol = 10^-p,cp = MachinePrecision},
While[(count++)<maxcount,
Block[{$MinPrecision = cp,$MaxPrecision=cp},
fx = f /. x->r;
dfx=df /. x->r;
step = fx/dfx;
r = r - step;
];
cp = -2*RealExponent[step];
If[cp ≤ MachinePrecision, cp = MachinePrecision,r = SetPrecision[r,cp]];
print[{count,step,cp}];
If[Abs[step] < Abs[r] tol,Break[]];
];
r
]

myRoot3[f1,{x,x01},30]

% - Sqrt[2]

Does the adaptive precision save much time?

print=.;
maxcount=100;

TableForm[data = Table[{p =2^k,First[AbsoluteTiming[r2 = myRoot2[f1,{x,x01}, p]]],First[AbsoluteTiming[r3 = myRoot3[f1,{x,x01}, p]]],RealExponent[r2-r3]},{k,4,20}]]

ListPlot[RealExponent[{data[[All,{1,2}]],data[[All,{1,3}]]}]]

Strength Training: Numerical Integration

Numerically compute (∫_0^100)sin((x^5))dx and (∫_0^100)sin((x^4) (Floor[x]))dx

NIntegrate[Sin[x^5],{x,0,100}]

NIntegrate[Sin[x^5],{x,0,100},Method->{Automatic,"SymbolicProcessing"->None}]

NIntegrate[Sin[x^4Floor[x]],{x,0,100}]

Block[{$MaxPiecewiseCases=200},NIntegrate[Sin[x^4Floor[x]],{x,0,100}]]

Levin Rule

Levin Rule applies to a large class of one-dimensional and multidimensional oscillatory integrals

Turns an oscillatory integral into non-oscillatory differential equations.

Applicable to:

Grid[{{"ExpRelated",Style[{Exp, Power}, "Input"]},
{"TrigRelated",Style[{Sin,Cos,Sinh,Cosh,Haversine,Sinc}, "Input"]},
{"ErfRelated" , Style[{Erf,Erfc,Erfi},"Input"]},
{"TrigIntegrals" ,Style[{FresnelS,FresnelC,SinIntegral,CosIntegral,SinhIntegral,CoshIntegral},"Input"]},{"BesselRelated",Style[{BesselJ,BesselY,HankelH1,BesselI,BesselK,HankelH2,StruveH,StruveL,SphericalBesselJ}, "Input"]},{"AiryRelated",Style[{AiryAi,AiryAiPrime,AiryBi,AiryBiPrime}, "Input"]},
{"HypergeometricRelated",Style[{HypergeometricPFQ,HypergeometricPFQRegularized,Hypergeometric0F1,Hypergeometric0F1Regularized,Hypergeometric2F1,Hypergeometric2F1Regularized},"Input"]},
{"InverseTrig", Style[{ArcSin,ArcCos,ArcSinh,ArcCosh,ArcSec,ArcCsc,ArcSech,ArcCsch,ArcTan,ArcCot,ArcTanh,ArcCoth,InverseHaversine},"Input"]},
{"Elliptic", Style[{EllipticE,EllipticK},"Input"]},
{"GammaRelated", Style[{Gamma,GammaRegularized,Beta,BetaRegularized},"Input"]},
{"ExpIntegrals",Style[{ExpIntegralE,ExpIntegralEi,DawsonF},"Input"]},
{"DifferentialRoots", Style[DifferentialRoot,"Input"]}, 
{"Other", Style[Log,"Input"]}}, Alignment->{Left,Center},Dividers->All]

Levin' s method applies to any sums or products of the above oscillatory functions, as well as their products and compositions with non - oscillatory functions.

NIntegrate[Sin[x^5],{x,0,100},Method->"LevinRule"]

Strength Training: Differential Algebraic Equations

Find the sensitivity of the solution to the following system of DAEs near a=2,b=-1/2

Clear[a];
daes={x'[t]+y[t]==Sin[2*Pi*t],y'[t]+z[t]==Cos[2*Pi*t],5*(a-Tanh[50*(t-1/2)])*y[t]+x[t]==-Sin[Pi*t],x[0]==-1/2};
dvars = {x[t],y[t],z[t]};
time = {t,0,1};

pf1=ParametricNDSolveValue[daes,dvars,time,{a}]

pf1[2]

pf2=ParametricNDSolveValue[daes,dvars,time,{a},Method->{"IndexReduction"->Automatic}]

pf2[2]

Plot[Evaluate[pf2[2]],{t,0,1},PlotRange->All]

Plot[Evaluate[pf2'[2]],{t,0,1},PlotRange->All]

Plot[Evaluate[(pf2[2.1]-pf2[1.9])/.2],{t,0,1},PlotRange->All]

Power Lifting: PDEs

Find the solution for the wave equation for a surface stretched over a disk where the initial shape is a cone

Using NDSolve directly

usol = NDSolveValue[{D[u[t,x,y],t,t]==Laplacian[u[t,x,y],{x,y}],u[0,x,y]==1-Sqrt[x^2+y^2],Derivative[1,0,0][u][0,x,y]==0,DirichletCondition[u[t,x,y]==0,True]},u,{t,0,5},{x,y}∈Disk[]]

plots = Table[Plot3D[usol[t,x,y],{x,y}∈Disk[],PlotRange->{-1,1},Mesh->None],{t,0,5,1/4}];

ListAnimate[plots]

Trying to be a bit more clever

The solution has to be axially symmetric, so the effective dimensionality can be reduced:

lu = Laplacian[u[t,r],{r,θ},"Polar"]

ϵ=0.1;
ursol = NDSolveValue[{D[u[t,r],t,t]==lu,u[0,r]==1-r,Derivative[1,0][u][0,r]==0,Derivative[0,1][u][t,ϵ]==0,u[t,1]==0},u,{t,0,5},{r,ϵ,1}]

ϵ=0.1;
ursol = NDSolveValue[{D[u[t,r],t,t]==lu,u[0,r]==If[r≤ 2 ϵ,1-2ϵ,1-r],Derivative[1,0][u][0,r]==0,Derivative[0,1][u][t,ϵ]==0,u[t,1]==0},u,{t,0,5},{r,ϵ,1}]

ϵ=0.1;
np = 100;
ursol = NDSolveValue[{D[u[t,r],t,t]==lu,u[0,r]==If[r≤ 2 ϵ,1-2ϵ,1-r],Derivative[1,0][u][0,r]==0,Derivative[0,1][u][t,ϵ]==0,u[t,1]==0},u,{t,0,5},{r,ϵ,1},Method->{"PDEDiscretization"->{"MethodOfLines","SpatialDiscretization"->{"TensorProductGrid", "MinPoints"->np,"MaxPoints"->np}}}]

rplots = Table[Plot[{usol[t,r,0],ursol[t,r]},{r,ϵ,1},PlotRange->{-1,1}],{t,0,5,1/4}];

ListAnimate[rplots]

Another way is to use the finite element for the spatilal discretization:

ϵ=0.001;
ursol = NDSolveValue[{D[u[t,r],t,t]==Derivative[0,2][u][t,r]+Derivative[0,1][u][t,r]/r,u[0,r]==1-r,Derivative[1,0][u][0,r]==0,u[t,1]==0},u,{t,0,5},{r,ϵ,1}, Method->{"PDEDiscretization"->{"MethodOfLines","SpatialDiscretization"->{"FiniteElement","MeshOptions"->"MaxCellMeasure"->Automatic}}}]

rplots = Table[Plot[{usol[t,r,0],ursol[t,r]},{r,ϵ,1},PlotRange->{-1,1}],{t,0,5,1/4}];

ListAnimate[rplots]

Summary

Hopefully, this workshop has given you some experience in working with numerical computation in Mathematica.

The range of things that can be done is vast, but once you understand a few principles, then it is easier to expand into new areas.

I strongly recommend looking at the documentation.  Not only are there lots of examples, but for some functions, particularly NDSolve, there are extensive tutorials explaining how the methods work and how best to use them with yet more examples: